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Abstract 

A general dispersion relation solver is presented which allows to numerical calculate all the solutions exactly and give 
polarizations naturally for magnetized fluid plasma with arbitrary number of components and arbitrary orient wave 
vector, with and without anisotropic pressure (e.g., firehose and mirror modes), relativistic, beam and gradient effects. 

Program summary 

Title of program: PDRF 
Catalogue identifier: 
Program summary URL: 

Program obtainable from: CPC Program Library, Queen University of Belfast, N. Ireland 

Computer for which the program is designed and others on which it has been tested: Computers: Any computer 
running MATLAB 7. Tested on DELL OptiPlex 380. 

Installations: Institute for Fusion Theory and Simulation, Zhejiang University, Hangzhou, PRC 

Operating systems under which the program has been tested: Windows XP Pro 

Programming language used: MATLAB 7 

Memory required to execute with typical data: 10 M 

No. of lines in distributed program, including test data, etc.: 340 

No. of bytes in distributed program, including test data, etc.: 12 000 

Distribution format: .tar.gz 

Nature of physical problem: Calculating all the solutions exactly and giving polarizations naturally for magnetized 
fluid plasma with arbitrary number of components and arbitrary orient wave vector, with and without anisotropic 
pressure, relativistic, beam and gradient effects. 
Method of solution: Solving as matrix eigenvalue problem. 

Restrictions on the complexity of the problem: Only linear multi-fluid, with temperature effects but without kinetic 
distribution function effects. 

Typical running time: About 1 second on a Intel Pentium 2.60 GHz PC. 
Unusual features of the program: Giving all the solutions and polarizations. 

Keywords: Plasma physics, Dispersion relation, Multi-Fluid, Waves, Instabilities, Matrix eigenvalue 
PACS: 52.21. Cm, 52.27 .Ny, 52.35.Qz, 52.35.Hr, 52.35.-g 



1. Introduction 

Since only few simple dispersion relations are analytical tractable in plasma physics, it is a historical issue to 
develop general numerical solvers for practical applications. Usually, we get the dispersion relation solutions from the 
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determinant of a 3-by-3 dielectric tensor, e.g, kinetic WHAMP (Waves in Homogeneous Anisotropic Multicomponent 
Magnetized Plasma) code by RonnmarkJTJ |2] and a fluid Mathematica code for magnetized parallel beam-plasma by 
BretG). 

It is difficult to generalize that usual treatment to arbitrary number of components with good convergence or to give 
all solutions. Here, we use a full-matrix approach to transform the task to an equivalent matrix eigenvalue problem. 
A general dispersion relation solver (PDRF, Plasma Dispersion Relation - Fluid version) for magnetized multi-fluid 
plasma with and without anisotropic, relativistic, beam and gradient effects, is developed. 

For space, astrophysical and laser plasma, this solver can have wide potential applications. 



2. Theory and Method 



We solve magnetized multi-fluid plasma, with Bq = (0, 0, Bq), Vjq - (Vjo x , Vjoy, Vjo z ), k = (k x , 0, k z ) - (k sin 9, 0, k cos 6) 
and Vrtjo/rtjo - (e„j X , e„j V , 0). We do not need derive the final 3-by-3 dispersion relation matrix for E = (E x , E v , E z ) as 
the usual treatment, such as by Stix[4 | and used by RonnmarklUEl and BretE). An alterant method is used: using the 
original full dispersion relation matrix and then treating the task as a matrix eigenvalue problem, instead of directly 
calculating its determinant. 

We consider only density gradient effects here and ignore the temperature gradient effects. 

2.1. Original equations 

Many (special and general) relativistic plasma models are just extensions of non-relativistic versions, and haven't 
been fully verified by experiments or observations. Some basic relativistic fluid plasma models can be found in |5|. 
We use usual choice 



d,rij = -V • (rijVj), 



VP, 



d,uj = -vj ■ Vuj + %(E + vjXB)-^f- Zfc - Vj)v ih 



(la) 
(lb) 
(lc) 
(Id) 



with u j 



d,E = c 2 V x B - J/cq, 
d,B = -V X E, 

'j ~ JjVj and 

J = XjqjtijVj, 

d,(P\\jpy v ) = o, 

'/■</'. //'/' ) = 0, 

which bases on [6| and [31, and where pj = mjnj, c 2 = l/fi^eo, jj = (1 - v 2 ./c 2 y 1 ^ 2 , and and y X j are the parallel 
and perpendicular adiabatic coefficients respectively. The treatment of anisotropic here is different from the CGL[7| 
one, but can reduce to the one by Bret and Deutsch[6| when setting = y ± j = y T j, where Py ± = nT\\ iX , with 
P = P\\bb + Pj_(I - bb) and b = B/B. When further setting T x ; - Ty, we get the isotropic pressure. 
After linearized, (|2]) gives 

J = IjQMjoVji+njiVjo), 



(2a) 
(2b) 
(2c) 



(3a) 



where c 2 = y ]Uj P ]Uj o/pjo and P j{) 



One should also note that 
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V ■ P fl = (ik x , 0, ik.) 



= c Lj n J 
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A jB xl ' 









AjB yl 


(4) 




AjB yl 


P\\fi . 





with Aj = (P\\jo - P± jo)/Bo and B^y j = 2^P\i X j I B 2 y where the off-diagonal terms come from tensor rotation (see 
similar expressions at Appendix A i from bo to b, are related to energy exchange and are important for anisotropic 
instabilities. A wrongly treating or ignoring of these off-diagonal terms will lose firehose or other important unstable 
anisotropic modes. 
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2.2. Full-matrix treatment 

The linearized version of Q with / = fo+f l e ihr - ib " is equivalent to a matrix eigenvalue problem (similar treatment 
can be found in |8| for MHD equations and [9| for ten-moment equations) 



AAX = MX, 



(5) 



with A - -icj the eigenvalue and corresponding X eigen vector, which represents the polarization information of each 
normal/eigen mode solution. 

For jlj, X = (n j uVju,Vjiy,Vju,E u ,Eiy,E u ,Bu,B 1 y,Bi z ) T , with« ; i = yp[v n +y 2 j0 (VjO • v n )v j() lc 2 ] = {a jpq }-v n 
(p,q = x,y,z) and y j0 = (1 - v 2 Q /c 2 )~ 1/2 . Here, A is 
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and, M is (vy terms when i + j haven't been written explicitly here) 
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y%v jQx Vjo y /c 
y%vj 0x vj 0z /c 2 



yjoVjOxVjOy/c 2 
y 3 j0 v j0zVj0y/c 2 



y%v j0x VjQ Z /c 2 

y%Vj0yVjOz/c 2 



yjo + y)o v )J c 



(1) 



(8) 



For s kinds of species, the dimensions of A and M are (4s + 6) x (4s + 6). Here, we can get all the solutions 
of the above system exactly (without convergence problem) via standard matrix eigenvalue solver, e.g. LAPACK or 
function eig() in MATLAB. 

MATLAB code PDRF is provided to solve the above eigenvalue problem. When giving all yj to 1, i.e., A — I, 



I, both be unit matrices, PDRF reduces to non-relativistic version. 



2.3. Breaks of linearization procedures 

During the linearization procedures, we have canceled all equilibrium quantities by assuming their effects are 
small, i.e., Q Q = Zjqjnjo - 0» Jo = Y. -'I}»M V » ~ °» v jo ' Vn >o - and F j0 = qjEo + qj(vp x B Q ) - VPjo/njo - 
m j Yji( v iO - v jo) v ij — 0. One should note that, usually, this is not the case in practices. For example, the non-zero Jo 
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Table 1: Breaks of linearization procedures. 
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Table 2: Comparing the cold plasma solutions using matrix method and Swanson's polynomial! 1 0"| . 



lo m ±10.5152 


±10.0031 


±9.5158 


±(2.4020E-4-i8E-19) 


±(1.1330E-4-ilE-16) ±0 ±0 


lj s ±10.5152 


±10.0031 


±9.5158 


±2.4020E-4 


±1.1330E-4 



will bring an extra Bqj field. To determined Bqj, we need boundary conditions or other information, which means 
that the system is not homogenous any more. More worse, these non-zero equilibrium quantities cannot be canceled 
by coordinate transformation, i.e., treating the system in moving frame. Another example is, for non-parallel beams, 
the Lorentz force will change the beams to ring beams at perpendicular direction during one cyclotron period. So, at 
most cases, we'd better use only parallel beams. However, for gradient modes or Eo + 0, there can exist non-zero Vj_o 
drift velocity in the system. 

All possible breaks of linearization procedures are shown in Table [T] One should keep these in mind when use 
PDRF. 



3. Benchmarks 



3.1. Cold plasma 

Cold limit (P s q = 0), without beam (v s o = 0), s = e, i, numerical solutions of {5} (w M ) and Swanson's fifth order 
polvnomial irTOll (a> s ) are given in Table|2j with kc = 0.1, 8 = n/3, = 1836 and u> pe = I0oj ce . The results are 

consistent with each other exactly, except some small (< 10~ 15 ) numerical errors. 

Scanning of k and 9 for co is shown in FigJT] with nti/m e = 4 and a> pe = 2co c 



3.2. Firehose and mirror mode 

Two well-known low-frequency hydromagnetic anisotropic instabilities are firehose mode and mirror mode. Roughly, 
the dispersion relation for firehose mode (or, anisotropic shear Alfven mode) is ^ 1 +Zj ^ ±J 2 ^" J • ^ or k±, using 
reduced expression from bi-Maxwellian kinetic dispersion relation, the mirror mode threshold is 1 +£7 P±j(l — j ± ) < 
ED. 

A correction for the term Aj in matrix element M V]U> b xI to 2Aj is needed to give the same mirror mode threshold as 
kinetic bi-Maxwellian plasma prediction since that there is an extra factor 2 in bi-Maxwellian plasma expression] 12 1, 
i.e., 

Spx = 2 P± (l ~ —)SB ± , (9) 
P\\ 

when comparing with Q. 

We check the instability thresholds of these two modes, with = and = 8, to scan fi ±i . A non-relativistic 
run is shown in Figj2] with and without the factor 2 correction for mirror mode. 

We can find that the numerical results consistent with analytical predictions very well. For the without correction 
case, the mirror mode threshold should be 1 + £j ^(1 - j^) < 0, which is confirmed in Fig Qa). When we use 
general CGL anisotropic fluid model, there will also exist other correction terms from SB to 6P. 
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Figure 1: Scanning of k and 8 for u>, with m\lm e = 4 and u) pe = 2ut l 




Figure 2: Scanning of f} ± for firehose and mirror modes thresholds. Dash lines are analytical predictions. 
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Table 3: Comparing the relativistic beam plasma solutions using matrix method and Bret's tensor| 3 1. 



oj m 3.2736 3.2731 0.9934 0.3000 0.3000 0.2890-i0.1664 0.2890+i0.1664 

B = u B 3.2736 3.2731 0.9934 0.3000 0.3000 0.2890-i0.1664 0.2890+i0.1664 

Un-magnetized "a? 3 E^16 0.0000 -0.0300 -0.0300 -1.0313 -3.2732 -3.2736 

oj b -0.0300 -0.0300 -1.0313 -3.2732 -3.2736 

J® 3.2945 3.2693 1.3427 0.5168 0.0771 0.2999+i0.0034 0.2999-i0.0034 

B * (jj b 3.2945 3.2693 1.3427 0.5168 0.0771 0.2999+i0.0034 0.2999-i0.0034 

Magnetized 0.0440 EA6 EA6 -0.1019 -1.3983 -3.2732 -3.2910 

oj b 0.0440 -0.1019 -1.3983 -3.2732 -3.2910 




Figure 3: Scanning of (k x , k z ) for relativistic electron beam mode with and without background Bq. 



3.3. Relativistic beam modes 

We test the default result of Bret's tensor[3| calculations. Ions immobile, or, no ions. Beam electrons y\ - 4.0, 
rib = O.lrtp. Background electrons v p = -vt,ni,ln p . For (Z v , 0,Z,) = kvb/a> pp = (0.3, 0, 3.0), PDRF matrix results a> M 
and Bret's tensor results a> B are listed in Table[3j which are exactly the same for both magnetized (Bo + 0, co ce - to pp ) 
and un-magnetized (Bo = 0) cases. 

Fig(3] shows the maximum growth rate beam mode at above parameters. To produce the same figure as Figj3] 
Bret's tensor in Mathematica takes about 1 minute, while PDRF takes a few seconds. 

3.4. The Doppler effects 

One should note that in the widely used non-relativistic multi-fluid (or kinetic) plasma equations, the field part is 
Lorentz invariant, while the fluid (or distribution function) part is Galilean invariant, then the total system is neither 
Galilean 

k' = k, a/ = a> — k ■ Vo, (10) 

nor Lorentz 

k' — y(k\\ - voco/c 2 ) + k ± , oj' — y(w - k ■ vo) (11) 

invariant, which will bring strange or incorrect Doppler effects, e.g., the solutions (both real frequency and growth 
rate) of electron beam and reverse ion beam would be very different. However, the relativistic version multi-fluid 
plasma equations should be Lorentz invariant. 

We check the Lorentz invariant ([TT} for two species (s = e, i) relativistic cold plasma mode. One should also 
note the relativistic density change no = y(vo)«Q, where vq is the moving velocity of the frame of reference. Initial 
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(a) co-k, 6=90°, cold wave LHW branch (b) |kx E|/|k- E|, n LH ^ ce n,.) 1/2 =0.5 




10 20 30 40 50 10 20 30 40 50 

kc kc 

Figure 4: Check the polarization of quasi-electrostatic LHW. 



parameters: q-, = —q e = 1, ni = n e = 4.0, Vo = 0, v,o z = Vo, v e o z — 0, k x c - 0.5, k z c = 1.0, which will give several 
solutions. We choose the largest real part solution a> = 2.792204183976196, and move to the ion frame of reference, 
i.e., v' iQy = but v^Q = -vo, which gives k' x c = 0.5, k'x - -3.471022809144551. The densities of electron and ion 
change~to n' eQ = 9.176629354822472 and n' iQ = 1.743559577416269. Using these new parameters as inputs, PDRF 
gives the largest real part solution u' M = 4.3410141 14998466, while (fTTli gives cj' l = 4.3410141 14998465. 

We can see that the only difference between a>' M and u' L is the last digit, which should be just numerical truncation 
error. 

For 3a> + mode, the Doppler shift wave vector is not limited to real number but is a complex number, which is 
also supported by PDRF. One can also check the effects of non-parallel beam, vy and P, to Doppler shift. However, it 
would be more complicated to verify them due to that more extra parameters transformations are needed. 



3.5. Low hybrid wave polarization and gradient drift instability 

Low hybrid wave (LHW, k\\ 0) is quasi-electrostatic mode, i.e., \kxE\ \ <§c \k-E\\. The polarization information 
of LHW solved by PDRF is shown in Figj4j which confirms the quasi-electrostatic conclusion. 

A further interesting benchmark is the low hybrid drift instability (reduced from [13] and Q/4|, Poisson equation 
is used) 

i -?-4^ + % i+ ^-^)- ' (i2) 

(a> — /cvo) (x> ce v k L c A ku> ' 

where e„ is gradient parameter and vo = v,o - v e rj is the drift velocity in electron frame of reference, which can be 
caused by density (or pressure) gradient or Eq x Bo. Using PDRF, one should firstly calculate the drift velocity, e.g., 
Vjo - VDj - _ 5~a;^( mn /o) x bo, as equilibrium quantity, for initial input. 



With Vq = and e„ = 0, ( 12 1 gives just the usual LHW solution, which is very close to PDRF solution, as shown 



in panels (a) and (b) in Fig|5]Setting e„, v = e ney — e„ — -0.15, v,o A = -v e Q X = vq/2 = 0.01, a typical result is shown in 
panels (c) and (d) in Fig|5] We can see that PDRF solutions also agree with ( fT2} qualitatively, but not quantitatively. 

This disagreement may be explained partly as: Dispersion relation ( 12 1 is just an approximate formula and which 
is also derived by Galilean transformation, while PDRF contains all effects and is not Galilean invariant. For example, 
when cold plasma with vq = but e„ + 0, the first row of matrix M is zero, which means that e„ won't affect the 
solutions at this case; while e„ will affect ( |T~2] > when vq = 0. 



3.6. More 

In this PDRF vl.0, the effects of pressure P, density gradient V«o/«o and collision v, ; , especially the relativistic 
effects of them, haven't been verified completely at present, mainly due to that it is difficult to find similar previous 
analytical or numerical results in literatures for comparisons. 
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However, it is easily for one to modify the model equations ([T} and |2]), or corresponding matrices A and M, if the 
present model breaks. 

4. Summary 

A general multi-fluid dispersion relation solver is provided, which shows good performances in several benchmark 
cases. This solver has the ability to calculate many more complicated effects, which are difficult or can not be solved 
by conventional solvers. 

The full-matrix treatment is a general and accuracy method for arbitrarily complicated fluid systems. However, 
kinetic effects will play important roles in many practical systems. While, the general treatment for kinetic dispersion 
relation with good convergence is still a difficult task in plasma community, especially for magnetized plasma, which 
contains n — oo orders cyclotron modes. 
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Appendix A. Extensions for general un-magnetized version 

Set Bo — in magnetized version of PDRF, we can get a simple un-magnetized version, e.g., as shown in Fig(3] 
In more general un-magnetized plasma, we need new treatment for anisotropic pressure relation, due to that there is 
no parallel background Bo- Here, for pressure and P x j, the parallel and perpendicular are to Vj when v ,o + 0. 
Without loss of generality, we can set k = (0, 0, k z ). 

A rotation matrix 
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is needed to transform each || ; and lj to x, y, z, i.e., the pressure tensor should be 



gives 
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(A3) 



For species Vjo = 0, one should set Pj = Pj as a scalar quantity. 

Other treatments are similar as the magnetized version. The above provides a guider for one who wants to develope 
a more general un-magnetized version of PDRF. 



Appendix B. FPDR User Manual 

PDRF contains two files: main program "pdrf.m" and input data file "pdrf.in". The input file has the follow 
structure 

qs ms ns vsx vsy vsz csz csp epsnjx epsnjy 

-1.0 1.0 4.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 

1.0 4.0 4.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 

One can add more species in new lines directly. 

Normalizations and other parameters are set in "pdrf.m". Since the code is not too long for understanding, one 
can modify default setup for new applications easily. Rewriting PDRF to FORTRAN or other languages is also very 
easy. 
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